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Abstract 

We are interested in an anisotropic singular diffusion equation in the plane and in its 
regularization. We establish existence, uniqueness and basic regularity of solutions to both 
equations. We construct explicit solutions showing the creation of facets, i.e. flat regions 
of solutions. By using the formula for solutions, we rigorously prove that both equations 
create ruled surfaces out of convex initial conditions as well as do not admit point (local) 
extrema. We present numerical experiments suggesting that the two flows seem not differ 
much. Possible applications to image reconstruction is pointed out, too. 



1 Introduction 

We study two examples of singular diffusion equations. One of them is an anisotropic total vari- 
ation (TV) flow, the other one is the same equation with the additive isotropic linear diffusion, 

* /Jdfrf^.^Y (1.1) 

| = 0^. + /3div(^,^). (1.2) 
at VKil \ u xi\J 

These problems are considered on a domain in R 2 . Our goal is to study features of solutions 
like facets, i.e. flat parts of solutions with normals corresponding to the singular directions. 
Our study was inspired by the phase transition theory appearing in crystal growth problems and 
image restoration, where presence of walls and edges plays a significant role, [13], [23], [24]. 

Let us describe ideas behind this note. The key element of the systems we study is the 
anisotropy. In both cases this determines the features of solutions. We will see that numerical 
experiments appear to give almost the same despite fact that the second equation is not degener- 
ate. The most spectacular phenomenon which is observed for this type of problems are flat parts 
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of solutions, connected with very strong diffusion, where Vw = 0. Such effects have been well 
studied for the isotropic total variation flow. We note that the interest in the TV flow arose from 
its application to image analysis and reconstruction, [25], [5], [2]. Namely, any regular level 
sets of solutions to this flow evolve by the mean curvature. This property is used to smooth out 
contours and in deblurring. We stress that numerical algorithms exploit properties of this flow 
even implicitly. 

The case of anisotropic diffusion is not so well studied. Despite the available papers like 
[18], [3], the mathematical theory is still far from the excellence. This changes however, because 
of the interest in algorithms detecting or retaining special image features like edges and corners. 
A conspicuous example is the paper [10] on 2D bar codes. We observe a growing body of 
literature devoted to this subject, [22], [7], [14], [16], [15]. We see the need to study evolution 
equations which are likely to preserve pronounced features of solution or its graphs besides 
facets, e.g. edges or corners. It turns out that the equations we study here may serve this 
purpose. The rigorous goals will be stated below, but we also present numerical simulations in 
section 5, which illustrate the qualitative features of solutions. 

We set the following main goals of this paper: 

£ to study facets, the flat parts of solutions, defined by Vw = 0; 

X to exhibit ruled surfaces, arising when one of the components of the gradient of solution u 
vanishes; 

& to construct special solutions, given by the explicit formulas which shows characteristic 
features of solutions; 

X to present numerical experiments and to show evolution of interesting model shapes, which 
were the motivation for looking for analytical results; 

9/t to propose a possible application to image processing. 

It is surprising that facets and ruled surfaces are the attributes of solutions to both systems. 
The lack of degeneration in the second model results only in smoothing out effect appearing 
near 'corners' and some dispersion. Hence in practice, one could find the second equation as 
more suitable for practical applications. Here, we present a series of solutions to both systems, 
represented by the gray scale. It shows some interesting differences, which nonetheless are very 
subtle. The initial data are represented by the last picture on Fig. 2. 

Applications to the phase transition theory are of particular interest, [24], however, systems 
(1.1) and (1.2) are rather a simplification of more complex models. In the image processing, 
usefulness may be more straightforward to see. The pictures presented in Section 5 show a 
possibility of reconstructing images. The diffusion in the second system helps us to restore the 
picture. Positivity of 7 gives averaging effects, but strong anisotropic nonlinearity keeps edges 
in the chosen directions. 

The upshot of these experiments is the following. The regularization of a very singular 
system yields not only smooth solutions but also it preserves the main features of the original 
equation. We should keep in mind this important observation in our future studies of systems 
with very singular nonlinear operators. 

At this point we note that the system with the added isotropic diffusion behaves like phase 
field models with respect to free boundary problems including the mean curvature flow. We 
mention just a few papers exploring the link, [1], [4], [6]. 
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Figure 1: Snapshots of evolution by (1.1) (first row) and (1.2) (second row) of the same initial 
datum 



Here, we do not plan to present a consistent theory, but rather to pinpoint a few interesting 
results and phenomena to find motivation for our future deeper mathematical analysis. In fact, 
this note can be viewed as an attempt to extend results for one-dimensional systems [17, 20, 
21, 19] on the two-dimensional case. To be more precise, we establish existence of solutions to 
both equations (1.1) and (1.2) by using the theory of nonlinear semigroups. For this purpose we 
exploit the gradient structure of (1.1) and (1.2). Uniqueness is automatically guaranteed. This is 
presented in the next section. There we also present exact formulas for solutions. The advantage 
is that they provide insight into the facet formation problem. Since the formulas do not always 
fit the framework of semigroup solutions, we recall the notion of a weak solution. The explicit 
solutions suggest that the flows of (1.1) and (1.2) make ruled surfaces out of the initial data, 
provided additional conditions are satisfied. This is rigorously established in Section 4. This 
considerations require quite precise regularity estimates established in Theorem in Section 3. 

The paper is organized as follows. In Section 2 we state basic existence results for systems 
(1.1) and (1.2) coming from the general theory. We point also to a few interesting explicit solu- 
tions illustrating typical shapes. In Section 3 we show that the solutions are of better regularity, 
provided the initial data are smooth enough. Next we prove conditional results, which explain 
why flat regions and ruled surfaces are typical for graphs of solutions. In Section 5 we con- 
centrate on numerical analysis and obtain a few interesting numerical solutions. These results 
show more direct phenomena which are able to be captured by the systems. In the appendix we 
present more complex example of an explicit solutions to (1.1). 
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2 Existence 



We will use general tools exploiting the structure of the problem. In order to use the semi- 
group theory we notice that we present equations (1.1), (1.2) as gradient flows of corresponding 
functionals on L 2 (f2). We set 

® o[U) ~ { +oo if u G L 2 (Q) \ BV(n). {lA) 

* (i a- J LilvuP + ^kJ + Ki) if^Gi/ 1 ^), 

$1 [U) ~ \ +oo if u e L 2 (fi) \ H 1 ^); (2 - 2) 

Here, f2 is an open subset of IR 2 , possibly unbounded, e.g. f2 = IR 2 . We study the above equa- 
tions on a square with homogeneous Neumann boundary conditions which is convenient from 
the numerical point of view, i.e. it is easier to implement a numerical scheme on rectangular 
domain. We also consider periodic boundary conditions. In general (3, 7 > 0, however we may 
scale the time and from now on we put (3 = 1 and admit 7 > or 7 = to have a possibility to 
study both cases simultaneously. 

It is obvious that $1 is well-defined and finite iff u G H l . The correctness of the definition 
of &o(u) is less obvious. In fact, this is an example of a more general situation studied in [18, 
Section 2]. Formula (2.1) should be understood as follows, 

/ (\u Xl \ + \u X2 \):=sup{[(z,Du)dx: z G Cq(Q; IR 2 ), ^ < 1}, (2.3) 
Jn Jn 

where |(pi,P2)|oo := max{|p!|, |_p» 2 1 }- It is now easy to check that these two functionals are 
convex, proper and lower semicontinuous on L 2 (VL). 

We notice that formally, the elliptic operator 7AM + div ^j^y , ^ j is the first variation of 

functional $1 while div (jf 3 ^, is the first variation of functional $ . Thus, equation (1.2) 
is the gradient flow of $1 and equation (1.1) is the gradient flow of $ . Keeping this in mind, 
we infer the following statement, where A^u) = —d&i(u), i = 0,1. 

Theorem 2.1 Let us suppose that u Q G D(Ai), i — 0,1. Then there exists a unique function 
u : [0, 00) — > L 2 {Vt) such that: 

(1) for all t > Owe have u(t) G D(A); 

(2 ) I e L-(0,oo,L 2 (fi)) and ||f \\ L ~ (0 ,oo,L*(n)) < K ? («o)IU^ 

(3) f G Ai(u(t)) a.e. on (0, 00); 

(4) u(0)=u . ' 

In addition, u has a right derivative at all t G [0, 00) and 

d + v 

+ A°(u(t)) = 0, (2.4) 

where A°(u(t)) is the minimal section of Ai(u{t)), (see [8]). 

Actually, since A,- L are subdifferential of convex functional, we say more. 
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Theorem 2.2 Let us suppose that u e L 2 and A4 are as in Theorem 2.1, then there exists a 
unique solution to equation 



d+u 



+ Ai( U (t)) 3 0, u(;0)=U . 



dt 



Moreover, for allt>0 u(t) belongs to D(A) and (2.4) holds. 

We notice that Theorem 2.1 follows from [8, Theorem 3.1], while [8, Theorem 3.2] implies 
our Theorem 2.2. In Theorem 2.1 we refer to the domains D(d& ) and D(d<&i), however, we 
abstain from exact description of theses sets. The semigroups obtained by these theorems are 
contraction semigroups, thus if Uq — > u in L 2 {Q), then for all fixed t we have u n (t) — > u(t). 
This observation will be used in the constructions of examples of solutions based on explicit 
calculations. 

This general result gives us the justification for our exact formulas for solutions. They are 
particularly valuable when we strive to study motion of facets or other special properties. First, 
for special data we cook up explicit formula for a solution to (1.1). To keep the simplest setting 
we consider the equations in the whole plane. We construct u (see formula (2.7)) a solution to 
a differential inclusion 



This initial condition does not belong to L 2 (IR 2 ), but u G L 2 oc (IR 2 ), u , Vw G BVi oc (R 2 ). The 
same property will be valid for u(- : t). Hence the notion of solution introduced in Theorem 2.1 
by (2.4) is not quite appropriate. This is why we introduce in (2.9) the notion of a weak solution. 

Proposition 2.1 Formula (2.7) below yields a weak solution to (1.1) in M. 2 with data (2.6), 
understood as (2.9). Moreover, (2.5) is satisfied in M. 2 x R + in a pointwise manner with the 
exception of a one dimensional set and the solution is Lipschitz continuous, but not C 1 . 

Proof. Let us define 



u t + A (u)3 ini 2 xl 



-+ 



(2.5) 



in place of (1.1), with the initial datum 



uq(xi, X2) = x\ + x\ — 2R 2 . 



(2.6) 



tit) ee t(t) = (I) 3 t l s ee -r (t) and hit) = (I) 3 tl+x 2 1+ x 2 2 - 2R 2 . 





We notice that for t > the quantities ^ (t) are uniquely defined by the condition 



o) = hit) = uoou^f)) with ea) = ±^h(t). 



The final observation is that these functions satisfy the equation 



dh 2 



h(0) = 0. 
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Now, we write the advertised formula for solutions to (1.1), 



u(x, t) 



2h(t) \xi\,\x2\<£(t), 
h(t)+x 2 -2R 2 \x 1 \<Z(t),\x 2 \>Z(t), 
h(t)+xj-2R 2 \x 2 \ < i(t), \xi\ > {(t), 



x 2 + x 2 - 2R 2 



(2.7) 



This formula defines a Lipschitz continuous, but not a C 1 function. 
We shall calculate u t . We obviously obtain 



u t (x,t) 



The point is to calculate a selection of C(Vu) := (u Xl /\u Xl \,u X2 /\u X2 \), where at least one 
of the components of Vw vanishes. For this purpose, we take advantage of the special structure 
of this operator, permitting us to use what we learned about the one dimensional case, see [20], 
[21]. This yields 



r 2h\t) 


X\ 


, k 2 | < €{t) 




\ h'(t) 


Xi 


<Z(t),\x 2 \ 




| h'{t) 


x 2 


<Z(t),\x 1 \ 


> m, 




Xi 


, W >£(t) 





C(Vu)(x,t) 



f^(x 1: x 2 ) |xi|, \x 2 \ < f(t), 

, sgnx 2 ) < £(t), \x 2 \ > £(t), 

sgnxi,|fy) \x2\<t(t),\x 1 \>£(t), 
'sgnxi, sgnx 2 ) \x 2 \ > £(t). 



( XI 



(2.8) 



This is a Lipschitz continuous vector field. Let us check if u is a weak solution to (1.1). We 
recall that u is a weak solution iff 



( Mt ,0) + (a,V0) = O inP'([0,T)) (2.9) 

for each e C™(R 2 x [0, T)) and a* e sgnu x ., i = 1, 2. 

Here we put a = £(Vtt)(i,t), where C(Vu)(x,t) is given by (2.8). Since a is Lipschitz 
continuous, we are allowed to integrate by parts in the second term of the LHS in (2.9), getting 
u t = div£(Vw). If we take into account the explicit form of h(t), it is easy to see that the 
identity holds everywhere, except a two-dimensional subset {(x 1: x 2 , t) : = £(£) or \x 2 \ — 
CW}ofR 2 xl + . □ 

This example was relatively easy to present, because the problem was consider on the whole 
M 2 . It is also interesting to see if a similar formula works on a bounded domain with a boundary 
condition. In Proposition 6.1 in the Appendix, we present a similar, but more messy formula 
for a square with Neumann boundary data. 

The same notion of weak solutions like introduced in (2.9) may be used also when Theorems 
2.1 and 2.2 are applicable. However, it is easy to see that if u satisfies (2.4), then it is a weak 
solution in the sense of (2.9). In addition, if u 1 and u 2 are weak solutions with the same initial 
data, then they must coincide. 

Next, we study solutions to (1.1) with data just in BV space. 
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Proposition 2.2 Let us suppose that Vt = (— L, L) 2 , and a G (0, L), M > 0. We set 



U (X 1 ,X 2 ) = -MX(-a,a)2. 

Then a unique solution to (1.1) with the above initial data is given by formula (2.10) below. 
Proof. Let us set 

otherwise (Z1U) 

Checking correctness requires defining £(Vtt) in a proper way. We define two auxiliary func- 
tions 

=^f xe(-L,-a), f -«zS^ xe(-L,-a), 

Z ~l( x ) = { \ x \ x \ ^ a > z *( x ) = { ~lT^ x M ^ a > 



fef a:e(a,L), I _Q! Z§=& xe(a,L), 



We now define £(Vu) by setting 

£(V«) = (Zi(xi)X{|x 2 |< a } + ^ 2 (a:i)X{|x 2 |>a}, ^l(a;2)X{|x 1 |<a} + ^(a^Xflzi^a})- 

It is now easy to check that 

u t = div£(Vw). 

We use the same argumentation as in the proof of Proposition 2.1, however the difference is 
that here the defined above C(Vu) is not so regular. In order to take the divergence we are 
required to control only appropriate directional derivatives, so the form of C(Vu) and Lipschitz 
continuity of Z 1 and Z 2 allow us to obtain the desired identity. This equality holds pointwise in 
M 2 x R+ except for a two-dimensional set. This formula is valid until the time when two facets 
merge into a constant stationary state at the extinction time t = T ext , 

2 2a 



T ext = M[- + — . □ 

Finally we point one special solution to the second system. We show existence of a moving 
front for (1.2), but without any boundary conditions. 

Proposition 2.3 Let us fix a > 0, then each of the functions given by the formula below is 
traveling front solution to (1.2), 



( 2t 

a 

2t , 1 J2 



u a (x,t) 



\xi\, \x 2 \ < a, 
+ ^x 2 < a, \x 2 \ > a, 

f + ^i \x 2 \ <a,\xi\> a, 

f + H^ + xl) \ Xl \,\x 2 \>a. 



Checking the correctness of this formula is easier than in the previous case. The above formula 
makes it clear that no traveling front solution is possible for (1.1). In the Appendix we point an 
extra explicit solution to (1.1). 
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3 Extra regularity 

In this part we show that solutions to (1.1) and (1.2) obtained via Theorems 2.1 and 2.2 are of a 
better regularity. It will be very important for deducing some qualitative features of solutions. 

Theorem 3.1 Let uo G H 1 ^), then the solution to (1.2) given by Theorem 2.1 fulfills the 
following estimate 

\\ut\\L 2 (o,TxK) + sup \\u t ^V 2 u\\ L ^ K) {t) < DATA{8). (3.1) 

te[s,T] 

Proof. We consider both cases at ones: 7 = and 7 > 0. After mollifying the system we 
test it by u t getting 

/ / ufdxdt + sup / [— |Vw| 2 + \u Xl \ + |u X2 |]dx 
Jo Jk te[o,T] Jk 2 

|«o,xil + \uo,x 2 \] dx - (3-2) 
The structure of the equation allows us to differentiate the system with respect to t. 

u tt - (d Xl (sgnu Xl ) t + d X2 (sgau X2 ) t + -fAu t ) = 0. (3.3) 

Let 77 be a time dependent function such that rj(0) = and for S > rj = 1, then we test (3.3) 
by u t r) getting 



sup / u 2 t r]dx+ / r][5(u Xl )u 2 Xlt + 5(u X2 )ul 2t + -f\Vu t \ 2 ]dxdt 
te[o,T] Jk Jo Jk 



< 2 f [ u 2 r]'dxdt. 
Jo Jk 

But the r.h.s. is bounded by (3.2), so we have u t G B(5,T; L 2 (K)). Taking into account the 
above information, we consider (1.2) in the following modification 

-[d xi (6(u xi )u xixi ) + d X2 (S(u X2 )u xlX2 ) + -fAu Xl ] = -u txi (3.4) 

here time is a fixed parameter. Testing (3.4) by u Xl , we get 

{S(u xl )u 2 XlXl +S(x 2 )u 2 XlX2 + >y\Vu xl \ 2 )dx < / \u t u xixi \dx 
k Jk 

which gives the estimates on 7 J K \ Vu Xl \ 2 dx. The same we have for x 2 . The estimate (3.1) is 
proved. □ 

If we use t 2 as a test function 77 above, then we obtain information on the blow up of ||« t ||. 
Namely, it is easy to see that 

Corollary 3.1 Under the assumptions of Theorem 3.1, we have ||mJl 2 < Ct~ x l 2 . 

We shall emphasize that the terms 5(u Xl )u 2 XlXl are considered just formal, to be precise we 
shall treat them as limits coming from analysis done on the level of approximation. 
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4 Ruled surface and convexity 



The first phenomenon, which is very expected for this type of systems, are features of minimiz- 
ers and maximizers of the solution. We ask about a possible structure of sets where the function 
u, for fixed time t, admits extrema. Since the issue of regularity is not well studied, here, we 
prove only the following result. 

Proposition 4.1 Let u be a solution to system (1.1) or ( 1.2). Let t > and for x in the domain 
u(-, t) has a minimum atx and in addition, u(-, t) is a convex function different from a constant 
in a neighborhood N of set u(-, t) = u(xq, t), then the set 

M = {x : u(x, t) = u(x , t)} fl N 

is a closed set with nonempty interior. 

Proof. We deduce that there is a sequence m n converging to m := u\ M from above and such 
that each level set {u(-, t) = m n } is a convex closed curve. Moreover, the sets M n = {«(•, t) < 
m n } are convex too. We integrate equation (1.2) over this set 



/ 



u t — ^Au — div( sgn-u^, sgn-u^J dx\dx2 = 0. 



Integration by parts leads us to the following conclusion, 

/ (77p + rii sgnu Xl + n 2 sgnn^dV, 1 = [ u t dx x dx 2 . 
J{ u =m n } on J Mn 

But convexity implies that |^ > at dM n . At the same time for almost all y functions x x H> 
u(xi,y, t) and x 2 u(y, x 2 ,t) are monotone, hence in a neighborhood of M n 

rii sgnw^ + n 2 sgnu X2 = \rii\ + \n 2 \ > \n\ = 1. 

We conclude that we obtain 

! dH 1 < \M n \ l / 2 ( ! uldx^) 1 ' 2 . 

J{u=m n } J M n 

Moreover, since u is not constant, then the sets M n must have a positive two-dimensional 
Lebesgue measure. However, due to the isoperimetric inequality we have 

U\dM n ) > ^—\M n \^ 2 , 

the identity holds for the disc. Hence 

C<(l u\dx x dx 2 ) 112 . (4.1) 

J M n 

However, due to Theorem 3.1, it* is square integrable, so the RHS of (4.1) above cannot go to 
zero when n — > oo. Thus, M is a convex set of positive two-dimensional measure, hence it 
must have nonempty interior. □ 
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The next feature concerns the shape of the graph of solutions. The example presented in the 
earlier section suggests that the graph of the solution develops parts which are ruled surfaces. 
To be more precise, we will show that if the level sets of a convex solution u(-,t) at t > are 
regular, then the graph contains ruled surfaces which are of positive two-dimensional measure. 
The tangent is orthogonal to vector (0, 1) or (1,0). The precise phenomenon is prescribed by 
the lemma below. 

Lemma 4.1 Let u be a sufficiently regular solution to (1.1) or ( 1.2), ( in other words 7 is equal 
to or 1). That means, for a fixed t the restriction ofu(-,t) to an open set U is convex. 
Furthermore, we assume that for given c£R, the level set 

S(c) = {x G K : u(t, x) = c} 

is regular, i.e. Vu\s( c ) exists T-l}-a.e. on S(c) and Vu|s( c ) 7^ "H 1 a.e. Then sets 





= {x: 


X 


Mf 


= {x : 


: x 


M+ 


= {x: 


X 


M 2 " 


= {x 


: x 



do not contain isolated points. 

Proof. It suffices to consider just one of these sets, e.g. M+. Let us suppose that our claim 
fails and M x + = {p}, in other words, function x 2 u(m + , x 2 ) has a strict minimum at x 2 = p 
in the interval [p — £,p + £}. Due to the continuity of u we notice that if u(xi, •) restricted to 
[p — £,p + £} attains its minimum on [p~(xi),p + (xi)], thenp~(x±), p + (xi) converge to p as x\ 
goes to m + . In particular, u(xi,p ± £) > u(m + ,p) for x\ close to m + . The last observation 
combined with monotonicity of u~ 2 (x 1: •), w+ (x 1: •) implies that 

u± 2 (x 1 ,p + e)>0, u± 2 { Xl ,p-£) < 0, 

for all xi close to m + . Thus, we can consistently define 



sgnu X2 



1 on {(x 1: p + £) : Xi G (m + — 5, m + + 5}, 
— 1 on {(^i,p — £) : Xi G (m + — 5, m + + 5} 



Let us take rectangles, R k = [m + , m + + 5k] x \p — £,p + £], where ^ < 5 and 5 is so small that 
the above considerations are valid. We integrate ( sgnu X2 ) X2 over R k . We obtain, 

/ (sgnu X2 ) X2 dxidx2 = / sgnM rr2 n 2 = 2 • 21 

We may assume that function x\ >->■ u(xi, p) is increasing on [m + , m + + <5], otherwise we could 
consider u(—x 1 ,x 2 ), in place of u(xi,x 2 ). 

Since x 2 h-> u(m,x 2 ) is convex, with minimum at a; 2 = p, then it must be increasing on 
[p,p + 5] and due to our assumption u(m + ,x 2 ) > u(m + ,p) for x 2 7^ p. Moreover, all lines 
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L = {(x, a) : x G M} intersect S(c) for a close to p, i.e. \p — a\ < 5, otherwise S(c) would be 
a point, i.e. a singular level set. Let us suppose 

(si,i 2 )el a n5(c), (4.2) 

with rr 2 close to p. Then, 

u±(m,x 2 ) > 0. 

Equality above is excluded because it contradicts (4.2) and the monotonicity of u xi (•, x 2 ). By 
the monotonicity of the derivative of a convex function we also obtain u xi (to, x 2 ) < u xi (to, x 2 + 
5 k ). Thus, we may consistently define sgnu Xl = 1 on the sides of R k parallel to the vertical 
axis. 

Let us now integrate our equation over R k , 

/ d Xl (sgau Xl ) + d X2 (sgnu X2 )dx 1 dx 2 = / (u t - -fAu)dx 1 dx 2 . 

jR k JR k 

performing integration by parts on the LHS and taking into account observations collected 
above, we conclude that 



/ sgn u X2 n 2 da = / (u t — ^An) dxidx 2 . 
JdRk JRk 



I dR k J R k 

We continue the calculations. Using the square integrability of u t —^Au established in Theorem 
3.1 we obtain that 

U < | f (u t - 7Aw) dx x dx 2 \ < (2£5 k ) 1/2 ( [ \u t - 7 Aw| 2 dx^\ 

jR k \JR k J 

i.e. 

A£ 1/2 < (25 k ) 1/2 Q \ut- 7 Am| 2 dx^j . 
If 5 k goes to zero, then we reach a contradiction. Thus, may not be a point. □ 

Theorem 4.1 Assume that for t > and a region A the solution u(-,t), restricted to A, is 
convex and the level sets ofu(-, t) satisfy the regularity assumption of Lemma 4.1, then sets 

Si = {(x, u(x, t)) : x e A, u Xl (x) = 0} and S 2 = {(x, u(x, t)) : x e A, u X2 (x) = 0} (4.3) 

are ruled surfaces, provided that m 4 ,7V 2 m is bounded pointwisely, and Vu is continuous for 
7 = 0. 

The proof of the above lemma follows immediately from Lemma 4.1. 
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5 Numerical experiments 



The algorithm used to perform numerical experiments is based on the duality approach consid- 
ered by Chambolle [9]. He computed a minimizer of the total variation model for the image 
denoising proposed by Rudin et al. [23]. In order to adapt this approach to solve the equation 
(1.2), we note first that the semi-discretization of (1.2) yields the following iterative scheme 



u m_ u m-l _ ^ _ ( u m „m 



= 7 A M m + /3V- , (5.1) 



where u m (x) := u(x,t m ) for m = 1,2,... and x £ M 2 , the initial data u°(x) := f(x) for 
x E IR 2 , where / e is a given function, and < St = t m — t m -\ for m — 1,2, ... 

denotes the time discretization step. For the convenience of notation, assume that St — 1 and 
consider the case m — 1. Then, we note that the equation (5.1) can be seen as the optimality 
condition for the minimization problem 



mm 



Q Jju- f) 2 + 1 \Vu\ 2 dx + (3 Jju xl \ + \u X2 \dx^j . (5.2) 



Let us introduce the differential operator A 1 : H l (Vi) — > defined by A^u := 

u — 7 Am. Using standard results of convex analysis (see, e.g., Ekeland and Temam [12]), we 
can show that the dual problem to (5.2) is 



mm 

gee, 



t$UlL A ^ f -^- 9Hf -^- 9)dx ) 



(5.3) 

subject to Igloo < 1, 



where g = (<?i, #2) is a vector function and |g|oo := max{|^i|, l^l}- 

From the Karush-Kuhn-Tucker conditions (see, e.g., Ciarlet [11, Theorem 9.2-4]), we get 
that there exist constants /ii, /j, 2 > 0, such that 

(A-\f-f3V-g)) Xk -n k g k = V, A; = 1,2, 

with either pL k > and \g k \ — 1 or = and \g k \ < 1 for k = 1,2. In any case, we have that 
A*i = \u xi I and /i 2 = \u X2 \, and therefore, we conclude that the solution u to problem (5.2) can 
be found by solving the system of equations 

, , ' (5.4) 

-«x fc + Pxjfi'ifc = , k = 1,2. 

In order to introduce the algorithm to solve (5.4), we need to turn into the discrete setting. 
From now on let f2 = (-L, L) 2 C IR 2 and values of the initial data / be given in the discrete 
set of iV 2 uniformly distributed points in fi. To simplify notation, we can fix the number N and 
take L such that N — 2L + 1. Now let / be a vector in the Euclidean space X = R N , defined 

by f(\x2 — L\ + l + \xi + L\N) := f(xi,x 2 ), for x±, x 2 — —L, —L + 1, L — 1, L, and let 
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us define vectors g±, g 2 and u in X in a similar way. Using this notation, we can introduce the 
discrete version of the system (5.4), given by 



A 7 u = f - (3 Yl=i D k9k , 

-D k u + \D k u\g k = , A; = 1,2, 



(5.5) 



where e Y with Y = ]R 7v2xAr2 i s a discrete version of the operator A 1 derived by the 
standard finite difference scheme taking into account the Neumann boundary conditions and 
(Di, D 2 ) G Y x Y corresponds to the discrete version of the gradient operator. To solve the last 
equations in (5.5), we follow Chambolle [9] and propose the fixed point iteration 

9k = 9t l + t (D k u n - \D k u n \g n k ) , k = 1, 2, 
for n — 1, 2, .... Finally, the algorithm to solve (5.4) is given by 



A,u- = f-/3 ELi D k g n k 



n-l 



9k 



9 k 



1 • v/V 



1 + r\D k u r - 



k = l,2, 



(5.6) 



for n — 1, 2, .... 



Theorem 5.1 Let r < (8Ai)~\ where \\ is the smallest eigenvalue of the operator A T Then, 
the sequence (u n , g n ) defined by the scheme (5.6) converges to the solution (u, g) of the equa- 
tions (5.5) as n — > oo. 

Proof. The proof can be carried out in a similar way as in Chambolle [9, Theorem 3.1] 
using the fact that A y is a symmetric positive define matrix, what implies that Ai > and 

(A^w, v) = (w, A~ l v ), for all w, v e X. 







Figure 2: Images f Sl , f Sa , f Ss and f Si . 



In the further part of this section, we present numerical solutions to the equations (1.1) 
and (1.2) with the Neumann boundary conditions and the initial data f s = — M%s, where 
Xs '■ (-L, L) 2 — » {0, 1} is a characteristic function of the set S C (— L, L) 2 . For experiments, 
we have taken L = 250, M = 50 and considered the following four sets: 

Si = {iGl 2 : ||x||i < 150} ,S 2 = {x ER 2 : \\x\\ 2 < 150} ,S 3 = {xeR 2 : Wx]]^ < 150} , 
S 4 = (Si U{xeR 2 : \\x - (125, 175) ||i < 25}) \{xGK 2 : \\x + (0, 125) ^ < 25} . 
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Images / 5l , fs 2 , fs 3 and f Si are presented in Figure 2. 

All experiments were performed with the same values for parameters involved in the algo- 
rithm, i.e., 7 = 5 _1 , /3 — 10, 8t — 1, r = 8 -1 . As the stopping criterion for the iterative scheme 
(5.6), we have used \\u n ^ 1 — w n ||2||w n ||2 1 < tol, with the tolerance tol = 1CT 5 . 



Figure 3: Numerical solutions to the equations (1.1) (upper row) and (1.2) (lower row) with the 
initial data f Sl , f Sa , f Sa and f Si , respectively. 



Numerical solutions to the equation (1.1) with the initial data accordingly equal to f Sl , fs 2 , 
fs 3 and fs 4 are presented in the upper row of Figure 3. The first two results have been obtained 
for m = 200 , whereas the next two results, for m = 170 and m = 90, respectively. We recall 
that m denotes the number of iteration of the scheme (5.1). Numerical solutions to the equation 
(1.2), with the same initial data and for the same numbers of iterations as before are presented 
in the lower row of Figure 3. 

The first two graphs in the upper row of Figure 4 present evolution of contour lines of 
solutions to the equation (1.1) with the initial data f Sl and f s , 2 , respectively. In each graph, 
contours are plotted for the level equal to the average value of a given initial data and correspond 
to solutions of the equation (1.1) for m — 0, 70, 140 and 210. The contour lines of solutions 
to the same equation but with the initial data f Ss and fs 4 and for m — 0, 60, 120 and 170 
are presented in the next two graphs in the same row. The lower row of Figure 4 presents the 
evolution of contour lines corresponding to solutions of the equation (1.2) with the same initial 
data and for the same numbers of iterations as before. 

The first two plots in Figure 5 show evolution of numerical solutions to the equations (1.1) 
and (1.2), respectively, along cross-sectional line x\ — passing through the middle of the 
square S\ for m — 0, 70, 140 and 210. In the case of the solution u to the equation (1.1), 
obtained values were equal to: {—5.25, —40.67} for m = 70, { — 10.5, —31.33} for m = 140, 
{—15.75, —22} for m = 210, where the first numbers in brackets correspond to values of u in 
Q\Si, and the second ones, in Si. We note that these results coincide with the exact values 
given by the formula (2.10) for t = f35tm. The third plot in Figure 5 presents the evolution 
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Figure 4: Evolution of contours corresponding to solutions to the equations (1.1) (upper row) 
and (1.2) (lower row) with the initial data f Sl , fs 2 , fs 3 and fs 4 , respectively. 




Figure 5: Evolution of numerical solutions to the equations (1.1) and (1.2) with the initial data 
fs 1 and to the equation (1.1) with the initial data G a * fs 1 . 



of the numerical solution to the equation (1.1) for smooth initial data obtained by convolution 
of the image fsj and the Gaussian kernel G a with the standard deviation a = 10. Similarly 
as in the one dimensional version of the equation (1.1) studied in [20]. Here, we may observe 
propagation of facets. 

In the last experiment, we were testing a possible application of the anisotropic total vari- 
ation flow equations (1.1) and (1.2) to solve the real problem of improving the quality of the 
scanned text. In this experiment, we were considering two binary images presented in the first 
column of Figure 6 with values scaled to {—50, 0}. Images in the second and third column of 
Figure 6 correspond to numerical solutions to the equations (1.1) and (1.2), respectively, for 
parameters 7 = 10~ 2 , f3 — 1, 8t — 1, r = 8~\ m = 15, and with images in the first column of 
Figure 6 as initial data. For comparison, in the last column of Figure 6, we present numerical 
solutions to the linear diffusion flow (the equation (1.2) for 7 = 10~ 2 , f3 — 0, St — 1, m — 15) 
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Figure 6: In columns: (a) initial data - two binary images of the scanned text, (b) solutions 
to the equation (1.1), (c) solutions to the equation (1.2), (d) solutions to the linear diffusion 
equation. 
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Figure 7: Results obtained by thresholding of images presented in the lower row of Figure 6 
on the level —10 (upper row) and on the level —5 (lower row). 



with the same initial data. We observe that in fact equation (1.2) represents the interplay be- 
tween an anisotropic total variation flow and the linear diffusion. It allows to fill corrupted parts 
of letters and at the same time slightly blur their boundaries. We notice that these properties are 
also visible in the results of the experiment with the image fs 4 , presented in the last columns of 
Figures 3 and 4. 
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In Figure 7, we present results obtained by thresholding images in the lower row of Figure 6 
on the level —10 and —5, respectively. We see that application of equation (1.1) gives basically 
better results, however in the case when larger parts of the letters are corrupted, the properties 
of equation (1.2) may be useful. In general, we infer from the experiments we performed that 
both total variation flow models analysed in this paper provide better results when applied to a 
class of real problem, than the standard linear diffusion equation. 



6 Appendix 

Formula (2.7) must be modified in order to accommodate the boundary conditions. This is done 
below. 



Proposition 6.1 Formula (6.1) below yields a weak solution to (1.1) in M. 2 with the data 

u (x 1 ,x 2 ) = {x\ + x\- 2R 2 )xb(o,r)(x 1 ,x 2 ) G L 2 (R 2 ). 

in the sense specified in Theorem 2.1. Moreover, the equation is satisfied in IR 2 in a pointwise 
manner with the exception of a one dimensional set and the solution is Lipschitz continuous, but 
not C 1 . 



Proof. Formula (2.7) shows the creation of a square facet and ruled surfaces over strips 
< £(0 and \x 2 \ < £(£). Now, we have to take into account their interaction with the 
boundary of the ball x\ + x\ < R 2 . The result is region Q(t), where u is different from zero. 
This set is defined as follows, Q(t) = B(0, R) n (-L(t), L(t)), where L(t) = ^/R 2 -f 2 (t). 

We shall see that the solution gets extinct, when the square facet hits the plane u = at 
t = ti. This is why for t G [0, ti), we set, 



' 2h(t) \xi 
h(t) + x 2 - 2R 2 \xi 
\xi 
u(x,t) = \ h{t)+x\ - 2R 2 \x 2 


x\ + x 2 - 2R 2 
I 



I, < £(0> O^i, £2) e Q(t), 
I < m <\x 2 \< y/R? -?(t), ( Xl ,x 2 ) e Q(t) 
I < at), N > y/B?- ?(t), {x u x 2 ) G Q(t), 

1 < at), at) <\xi\< y/n?-e(t), ( Xl , X2 ) & m, 

\<at),\xi\ > y/B?- ?(t), 

_|, \x 2 \ > at), (xi,x 2 ) G Q(t), 



x 2 

Xi 



This formula is valid up to 2£ 2 (ti) = R 2 , i.e. U = ^R 3 . 

Calculating Vm is easy, but we have to modify C(Vu). Namely, we set, 



(6.1) 



C(Vu) 2 = 





( XI 

m 


\ x i\ 


-1 


sgnxi 






I 


\x 2 \ 




f %2 

m 


\x 2 \ 


-\ 


sgnx 2 


\x 2 \ 




{ 


\ x l\ 



<at),\*i\ < y/R 2 -?(t), 
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We notice that vector field C(Vu) has jump discontinuities, nonetheless its distributional di- 
vergence is in Lf oc and has the desired properties. It is now easy to check that u satisfies (1.1) 
pointwise except a two-dimensional set in IR 2 x R + . We note the discontinuity of u t is respon- 
sible for the creation of the two dimensional facet and its growth. □ 
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